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ABSTRACT 

We investigate numerically and semi-analytically the collapse of low-mass, rotating 
prestellar cores. Initially, the cores are in approximate equilibrium with low rotation 
(the initial ratio of thermal to gravitational energy is ao — 0.5, and the initial ratio 
of rotational to gravitational energy is /3rj = 0.02 to 0.05). They are then subjected to 
a steady increase in external pressure. Fragmentation is promoted - in the sense that 
more protostars are formed - both by more rapid compression, and by higher rotation 
(larger f3 ). 

| In general, the large-scale collapse is non- homologous, and follows the pattern 

. described in Paper I for non-rotating clouds, viz. a compression wave is driven into 

the cloud, thereby increasing the density and the inflow velocity. The effects of rotation 
become important at the centre, where the material with low angular momentum forms 
a central primary protostar (CPP), whilst the material with higher angular momentum 
forms an accretion disc around the CPP. More rapid compression drives a stronger 
compression wave and delivers material more rapidly into the outer parts of the disc. 
Consequently, (i) there is more mass in the outer parts of the disc; (ii) the outer parts 
of the disc are denser (because the density of the material running into the accretion 
shock at the edge of the disc is higher); and (iii) there is less time for the gravitational 
torques associated with symmetry breaking to redistribute angular momentum and 
thereby facilitate accretion onto the CPP. The combination of a massive, dense outer 

• disc and a relatively low-mass CPP renders the disc unstable against fragmentation, 

and leads to the formation of one or more secondary protostars. At their inception, 
these secondary protostars are typically four or five times less massive than the CPP. 

For very rapid compression there is no CPP and the disc becomes more like a 
ring, which then fragments into two or three protostars of comparable mass. 

For more rapid rotation (larger (3q), the outer disc is even more massive in com- 
parison to the CPP, even more extended, and therefore even more prone to fragment. 
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1 INTRODUCTION 

It is now well established that stars form in dense cores 
embedded in molecular clouds (see e.g. Andre et al. 2000). 
However the conditions under which these cores form, be- 
come unstable, collapse and fragment are still a matter of 
debate. 

In a previous study, Hennebelle et al. (2003, Paper 
I) have investigated the possibility that the collapse of a 
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prestellar core is driven from the outside by an increase in 
the external pressure. This model seems to reproduce many 
of the key features observed in nearby star forming cores. 
In particular, (i) during the pre-stellar phase, the density 
profile is approximately flat in the centre (Abergel et al. 
1996, Bacmann et al. 2000, Motte & Andre 2001). (ii) For 
slow to moderate compression rates, subsonic infall veloci- 
ties develop in the outer parts of the core during the prestel- 
lar phase (Tafalla et al. 1998, Williams et al. 1999, Lee et 
al. 1999). (iii) During the Class phase, subsonic veloci- 
ties persist in the outer parts of the core (Belloche et al. 
2002), and transsonic velocities develop in the inner parts 
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(r<2000AU). (iv) There is an initial short phase of rapid 
accretion onto the central protostar (the Class phase), fol- 
lowed by a longer phase of slower accretion (the Class I 
phase), as inferred from the studies of Greene et al. (1994), 
Kenyon & Hartmann (1995), Bontemps et al. (1996) and 
Motte & Andre (2001). 

A question that was not addressed in Paper I is the 
formation of multiple protostars by fragmentation of a col- 
lapsing core. This question is critical, because it is now well 
established (Duquennoy & Mayor 1991, Fischer & Marcy 
1992, Chez et al. 1997) that a large fraction of stars is in 
binary - or higher multiple - systems (i.e. at least 50% for 
mature solar-type stars in the field, and a higher percentage 
for Pre-Main Sequence stars in young associations). 

1.1 Previous Work 

The gravitational fragmentation of a collapsing cloud has 
been intensively investigated for over two decades (see Bo- 
denheimer et al. 2000 for a review), mostly with numerical 
simulations, but also by semi-analytic means. The enduring 
hope has been to derive a robust general theorem defining 
the conditions required for fragmentation, for instance, in 
terms of the initial ratio of thermal to gravitational energy, 
Qo , and the initial ratio of rotational to gravitational energy, 
/3o- However, robust theorems have proved elusive. This is 
in part because the parameter space of initial conditions 
and constitutive physics is very large. Therefore the gener- 
ality of the results obtained is hard to establish. In addition, 
fragmentation appears to depend sensitively on the compli- 
cated radiative-transport and thermal-inertia effects which 
come into play when star-forming gas switches from being 
approximately isothermal to being approximately adiabatic. 
The computational resources required to simulate this radia- 
tive transport properly are not yet available. Furthermore, 
all such theorems inevitably beg the question: How were the 
unstable initial conditions created in the first place? Here 
one probably has to appeal loosely to the chaotic effects of 
supersonic turbulence in star forming molecular clouds (e.g. 
Elmegreen 2000; Padoan & Nordlund 2002) 

The simplest case to treat is a spherical cloud with uni- 
form density, solid-body rotation, and isothermal equation 
of state; this is sometimes called the standard model. Us- 
ing semi-analytic arguments, based on the classical stabil- 
ity analysis of Maclaurin spheroids (Lyttleton 1953; Chan- 
drasekhar 1969), Tohline (1981) concluded that all such 
clouds should fragment. Subsequently, Miyama, Hayashi 
& Narita (1984) used numerical simulations to derive a 
fragmentation criterion of the form ao/3o^0.12. A simi- 
lar criterion was obtained semi-analytically by Hachisu & 
Eriguchi (1984, 1985), and Miyama (1992) revised the crite- 
rion slightly to a?oA) J$ 0.15, on the basis of an analysis of the 
flatness and stability of rotating discs. Tsuribe & Inutsuka 
(1999a), again using a semi-analytic approach, took account 
of the non-homologous nature of collapse, and obtained a 
somewhat different criterion, which can be approximated 
by a a < 0.55 — 0.65/3o- This criterion was subsequently con- 
firmed with simulations (Tsuribe & Inutsuka 1999b). 

The standard model can be modified in several inter- 
esting ways. 

(1) The standard model can be modified by adding 
an m = 2 azimuthal density perturbation, i.e. p(r) = 



po [1 + Acos((f>)], where <j) is the azimuthal angle in spher- 
ical polar co-ordinates and A is the fractional amplitude of 
the perturbation. This perturbation is a common basis for 
numerical explorations of collapse and fragmentation. The 
standard model was first simulated by Boss & Bodenheimer 
(1979), who invoked a perturbation with amplitude A = 0.5. 
They found that a cloud having Qo = 0.25 and flo = 0.20 
collapsed and fragmented to form a binary system. This re- 
sult was reproduced by Burkert & Bodenheimer (1993) using 
a more accurate code. Burkert & Bodenheimer (1993) also 
performed simulations with ao = 0.26 and j3o = 0.16 and 
with the amplitude of the perturbation set to A = 0.5 and 
A = 0.1. In the latter case they obtained not just a binary, 
but also a line of smaller fragments between the main binary 
components. However, Truelove et al. (1998) have demon- 
strated using an AMR code, that this is an artefact. If the 
collapse remains truly isothermal to high densities, and if 
the resolution is sufficiently high for the Jeans mass to be 
resolved at all times, the material between the binary com- 
ponents should not fragment. This has been confirmed by 
Boss et al. (2000), by Sigalotti & Klapp (2001), and by Kit- 
sionas & Whitworth (2002), who followed the purely isother- 
mal collapse to even higher density, using SPH with particle 
splitting. 

(2) The standard model can be modified by changing 
the equation of state. Tohline (1981) and Miyama (1992) 
considered the collapse of clouds having an adiabatic equa- 
tion of state, P = Kp' 1 , and concluded that the condition 
for fragmentation takes the form 

«(4 4 - 37) < /(7) , (1) 

where, for example, /(5/3) ~ 0.064. Alternatively, a 
barotropic equation of state can be adopted, in which the 
gas is isothermal at low densities (where star forming gas 
is expected to be thin to its own cooling radiation), and 
adiabatic at high densities (where the gas is expected to 
be optically thick to its cooling radiation). Collapse simula- 
tions using a barotropic equation of state of this form are 
presented by Bonnell (1994), Bate & Burkert (1997), Boss 
et al. (2000), and Cha & Whitworth (2003). A barotropic 
equation of state of this type is also used in the simulations 
presented in this paper, and is discussed in Section 12.21 It 
is clear that treatment of the equation of state has a pro- 
found influence on the outcome of collapse. For example, 
Bate & Burkert (1997) show that the collapse of an initially 
uniform-density cloud having ao = 0.26, /3o = 0.16, and an 
m = 2, A = 0.1 azimuthal perturbation, does produce a 
line of small fragments between the two main binary com- 
ponents, if the barotropic equation of state described above 
is used (but not if the gas remains isothermal indefinitely). 
Perhaps the most telling result in this regard is reported by 
Boss et al. (2000) who simulated the collapse of a cloud with 
an m = 2, A = 0.1 perturbation, first using a barotropic 
equation of state, and then including radiation transport 
and an energy equation. The latter case produced a binary, 
whereas the former case did not, even though the variation 
of pressure with density was very similar in the two cases. 
The implication is that the subtle radiation-transport and 
thermal-inertia effects, which suddenly become important as 
isothermality gives way to adiabaticity, play a critical role 
in determining the pattern of fragmentation. 

(3) The standard model can be modified by imposing a 
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density profile. Myhill & Kaula (1993), using a code which 
included radiation transport and an energy equation, showed 
that clouds with solid-body rotation, ao = 0.16, /3o = 0.17, 
and an m — 2, A = 0.1 or 0.5 azimuthal perturbation do 
not fragment during the isothermal collapse phase, if the 
density profile is centrally peaked (specifically p cx r~ n with 
n = 1 or n = 2). However, Burkert, Bate & Bodenheimer 
(1997) repeated the simulation with A = 0.1 and n = 1, 
using a barotropic equation of state. They found that after 
the central primary protostar condenses out, a circumstellar 
disc forms around it, and this disc then fragments to produce 
companion protostars. A number of simulations have also 
been performed with Gaussian density profiles (e.g. Boss 
& Myhill 1995; Boss 1996; Burkert & Bodenheimer 1996; 
Boss et al. 2000), and also with exponential profiles (Boss 
1993), but there does not appear to have been a systematic 
evaluation of the influence of these profiles on the outcome 
of collapse. 

(4) The standard model can be modified by introduc- 
ing differential rotation. Myhill & Kaula (1993), again using 
a code which included radiation transport and an energy 
equation, showed that clouds with ao = 0.16, /3q = 0.17, 
a centrally peaked density profile (p oc r _n with n = 1 or 
n = 2), and an m — 2, A = 0.1 or 0.5 azimuthal perturba- 
tion, do fragment if they have sufficient differential rotation 
(in contrast with the result for solid-body rotation). The 
tendency for differential rotation to promote fragmentation 
has been confirmed by Boss & Myhill (1995) and Cha & 
Whitworth (2003). 

(5) Finally, a number of authors have explored the ef- 
fect of clouds having non-spherical initial shapes. Evidently 
the outcome here depends not only on ao (and f3o, if there is 
rotation), but also on the aspect ratio of the initial configu- 
ration. Bastien (1993) simulated the collapse of isothermal, 
non-rotating cylindrical clouds and determined the mass per 
unit length required for fragmentation. This problem was 
revisited by Bonnell & Bastien (1991) and Bastien et al. 
(1991), and extended to polytropic cylinders by Arcoragi 
et al. (1991). Bonnell et al. (1991) explored the circum- 
stances under which isothermal cylinders rotating about an 
axis perpendicular to their elongation fragment to produce 
binary systems, and Bonnell et al. (1992) considered isother- 
mal cylinders rotating about an arbitrary axis. Bonnell & 
Bastien (1992) repeated this last study with cylinders hav- 
ing a density gradient along the symmetry axis, and showed 
that quite modest gradients were sufficient to produce bi- 
nary systems with mass ratios in the range 0.1 to 1. Nel- 
son & Papaloizou (1993) simulated the collapse of prolate 
spheroids, and showed that they fragment if the mass per 
unit length is sufficiently high (as for cylinders), but the bi- 
nary components tend to be closer because there is less mass 
at the ends of a prolate spheroid. Boss (1993) simulated 
the collapse and fragmentation of a rotating, mildly prolate 
cloud with an exponential density profile, and showed that 
the conditions for fragmentation were rather more restrictive 
than those obtained by Miyama, Hayashi and Narita (1984). 
These results were extended to mildly prolate clouds having 
a Gaussian density profile and differential rotation by Boss 
& Myhill (1995); and to more elongated prolate clouds hav- 
ing a Gaussian density profile, but solid-body rotation, by 
Sigalotti & Klapp (1997). Boss (1996) simulated the collapse 
of rotating isothermal oblate spheroidal clouds and showed 



that fragmentation into multiple systems occurs provided 
Qo < 0.4, almost independent of /3q. Monaghan (1994) ex- 
plored the effect of vorticity on the collapse and fragmenta- 
tion of ellipsoidal clouds, using numerical simulations. 

1.2 Outline of paper 

In this paper we pursue further our investigation of prestel- 
lar cores whose collapse is induced by a steady increase in 
the external pressure. We extend the investigation to the 
case of rotating cores, and focus on the formation of circum- 
stellar discs around the central primary protostars (CPPs) 
and the subsequent fragmentation of these discs. Our main 
goal is to relate the properties of the star systems formed 
(multiplicities and mass ratios) to the dynamics of the col- 
lapse, and hence to the parameters f3o and <f>, measuring - 
respectively - the initial rate of rotation and the rate of com- 
pression (see Eq. @). Our model invokes initial conditions, 
and generates density and velocity profiles, very similar to 
those inferred from observations of real star forming cores. 
However, without performing a large (i.e. statistically ro- 
bust) ensemble of simulations, we cannot know whether it 
also delivers multiple systems with statistics similar to those 
observed. 

In Section [5] we describe the numerical method we use, 
the constitutive physics, and the initial and boundary con- 
ditions. In Section |3] we present our results, with special 
attention to the large-scale velocity and density fields in the 
cores. Section [I] discusses the evolution of the discs that 
form around the central primary protostars, with particular 
emphasis on the fragmentation process. Section summa- 
rizes our main conclusions. In the Appendix, we develop a 
semi-analytical description of the key features discussed in 
Sections [3] and 2] an d estimate the timescales influencing 
disc stability. We perceive this analysis as an integral part 
of our paper, and we have put it in an appendix purely so 
that those whose main interest is in the phenomenology of 
fragmentation can read about the results of our simulations 
without getting to grips with the more mathematical as- 
pects, which give an approximate quantitative explanation 
for the dependence of disc stability on the initial core rota- 
tion and the rate of compression. 

2 CONSTITUTIVE PHYSICS, INITIAL 

CONDITIONS AND NUMERICAL METHOD 

2.1 Co-ordinates 

With respect to a Cartesian co-ordinate system (x, y, z), the 
global angular momentum of the core will always be directed 
along the z-axis. Distance from the origin will be denoted 
by 



and distance from the rotation axis by 

w = [x 2 + y 2 ] 1/2 . (3) 

The velocity is then divided into three orthogonal compo- 
nents: the equatorial velocity component, 

v w = [xv x + yv y ]/w ; (4) 
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the azimuthal velocity component, 

Vg = [xv y - yv x ]/w; (5) 
and the polar velocity component, v z . 

2.2 The Equation of State 

We use a barotropic equation of state (cf. Bonnell 1994), 
which mimics the expected thermal behaviour of star form- 
ing gas (e.g. Tohline 1982, Masunaga & Inutsuka 2000): 
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Here P is the pressure, p is the density, and C s is the isother- 
mal sound speed. 

At low densities, p < po = 10~ 13 gcm~ 3 , C s — Co — 
0.19 kms - , corresponding to isothermal molecular gas at 
10 K. The presumption is that the gas is able to radiate 
freely, either via molecular line radiation, or - once the den- 
sity rises above about 10 -19 gem -3 - by coupling thermally 
to the dust. 

At high densities, p > po, P oc p 5 ^ 3 , corresponding to 
an adiabatic gas with adiabatic exponent 7 = 5/3. Here 
the presumption is that the cooling radiation is trapped by 
dust opacity. We note that molecular hydrogen behaves like 
a monatomic gas until the temperature reaches several hun- 
dred Kelvin, because the rotational degrees of freedom are 
not excited at lower temperatures, and hence 7 = 5/3 is the 
appropriate adiabatic exponent. 

The switch to adiabatic behaviour at high density 
causes the gravitational collapse to slow down, and obviates 
the need to invoke sink particles. This allows us to capture 
the dynamics of the disc, and accretion from the disc onto 
the central protostar, much more accurately. If a sink were 
introduced it might seriously influence the development of 
spiral density waves in the disc and its patterns of fragmen- 
tation. 



2.3 The Initial and Boundary Conditions 

The initial conditions are a rotating truncated Bonnor-Ebert 
sphere, contained by a hot rarefied intercore medium having 
uniform density. The Bonnor-Ebert sphere is truncated at 
£ = 6 (i.e. at radius R = 6Go/(47rGp c ) 1/2 , where G is the 
universal gravitational constant and p c is the central den- 
sity). The core mass is one solar mass and its initial radius 
is ~ 0.05 pc. This is a reasonable representation of observed 
prestellar cores (e.g. Andre, Ward- Thompson & Barsony, 
2000). 

Because of its initial rotation, this configuration is not 
strictly in equilibrium. However, since the rotational energy 
is only a few percent of the gravitational potential energy 
(/3q = 0.02 to 0.05), it is very close to equilibrium. 

Molecular-line observations of dense cores (Goodman 
et al. 1993) suggest that typically the ratio of rotational to 
gravitational energy is 0o — 0.02. We therefore adopt this 
value for most of our simulations. In addition, we consider 
/3o = 0.05, in order to explore the dependence on /3o- 

The rotation profile is obtained by assuming that the 
original core had uniform density and rotated as a solid- 
body, and that angular momentum was then conserved 



minutely whilst the core evolved from this original uniform- 
density state to our centrally condensed initial conditions. 
This means that our initial configuration is rotating differ- 
entially at the outset. 

At t — 0, the temperature of the intercore medium is 
increased in such a way that its pressure satisfies 



P{t) = Po + Pt. 
We define 
, _ Po/P 



(7) 



(8) 



r Ro/Co ' 

Thus 4> is the ratio between the time scale on which the 
intercore pressure doubles and the sound- crossing time of 
the core. A low value of <j> means rapid compression, whereas 
a high value means slow compression. 

2.4 Numerical method 

The numerical method used is very similar to that described 
in Paper I. We use a standard Smoothed Particle Hydro- 
dynamics code (e.g. Monaghan 1992) in combination with 
Tree-Code Gravity. There are three types of particle: the 
core particles, which experience both hydrodynamic and 
gravitational forces; the intercore particles, which experi- 
ence only hydrodynamic forces; and the boundary particles, 
which are passive. There are ~ 10 5 core particles, ~ 5 x 10 4 
intercore particles, and ~3x 10 4 boundary particles. 

Since we are simulating rotating cores, the intercore and 
boundary particles are given an initial uniform angular ve- 
locity, so as to minimize loss of angular momentum due to 
friction between the core and intercore gas, and between the 
intercore gas and the boundary. 

2.5 The Jeans condition 

In a numerical simulation involving self-gravity, it is essential 
that the Jeans condition be obeyed, i.e. that the Jeans mass 
be resolved at all times (Bate & Burkert 1997; Truelove et al. 
1997, 1998). The Jeans mass is Af. Jcans ~ 6G _3/2 p" 1/2 C* 3 , 
and the minimum resolvable mass in SPH is M reso ived ~ 
A/"neibfi, where A/" nc ib — 50 is the mean number of neigh- 
bours within the smoothing kernel of a typical SPH particle, 
and m is the mass of a single SPH particle. Thus the Jeans 
condition can be written as an upper limit on the mass of a 
single SPH particle, 

6G 3 

m<m max ~ A/ . eibG / /v/2 . (9) 

From Eqn. © we see that the minimum value of C 3 /p 1 ^ 2 is 
2 3//2 Cq/Pq 2 . Hence the Jeans condition is always satisfied 
as long as 

2 3//2 fiC 3 

m < m max : — ^ ~ 2.5 x 10~ 4 M Q . (10) 



1/2 



A/; cib G 3 / 2 Po 

Since we model a 1M core with 10 equal-mass particles, 
we have m = 10 _5 Mq. Consequently the Jeans condition 
is easily satisfied, and this ensures that the fragmentation 
which occurs in our simulations is not a consequence of poor 
numerical resolution. In addition, we have repeated all the 
simulations presented here with ~ 5 x 10 4 core particles, and 
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shown that the results are statistically unchanged; precise 
agreement is not expected, since fragmentation is seeded by 
particle noise, and particle noise is dependent on the number 
of particles. 



2.6 Angular momentum conservation 

Although the SPH equations ensure conservation of the glo- 
gal angular momentum, angular momentum is not necessar- 
ily well conserved locally. In particular, for the simulations 
involving slow compression, the duration of the pre-collapse 
phase can be as long as three or four freefall times, and sig- 
nificant non-physical transport of angular momentum can 
occur during this time. By non-physical transport of angu- 
lar momentum we mean the transport which arises before 
azimuthal symmetry is broken, due to differential rotation 
and the friction caused by artificial and numerical viscosity. 

In order to limit this effect, we invoke the Balsara switch 
(Balsara 1995), i.e. we multiply the artificial viscosity term 
by the factor 



IV.vl 



|V.v| + |V A v| +C 3 /(10 3 h) 



(11) 



Even with this factor fEquation lllL we find that angu- 
lar momentum loss from the inner parts of the core can be 
significant. For example, with the slowest compression rate 
that we treat (<f> = 3, see next section), the 10,000 densest 
particles (10% of the total number of particles) have lost 
about 20% of their initial angular momentum by the time 
the inner disk starts to form; in contrast, the total cloud 
angular momentum is conserved to within a few percent. 
For the intermediate compression rate (4> = 1), the 10,000 
densest particles have lost about 10% of their initial angular 
momentum by the time the inner disc starts to form, and 
for the faster compressions (0<O.3), this figure is <5%. 

In order to check that our results are not significantly al- 
tered by non-physical transport of angular momentum, we 
have repeated several simulations with local conservation 
of angular momentum imposed on all particles having den- 
sity p < 10 _4 po- This ensures that angular momentum is 
conserved locally as long as the core remains axisymmet- 
ric. Physical transport of angular momentum, due to the 
gravitational torques which accompany symmetry-breaking 
instabilities, does not occur until the disc density exceeds 
10 _4 po (see Section 3J. 

In order to impose local conservation of angular momen- 
tum, at each time-step and for each particle i, we calculate 
the velocity by solving the equation of motion, and then 
we extract the azimuthal component, v$,i{t). ve,i{t) is then 
recalculated so as to enforce local conservation of angular 
momemtum, i.e. 



ve,i(t) = v eii (0)wi(Q)/wi(t). 



(12) 



In these simulations the angular momentum loss of the 
10,000 densest particles is reduced to about 5%, before sym- 
metry breaking occurs. However, in terms of the growth and 
fragmentation of the central disc, there is no significant dif- 
ference from the simulations where local conservation of an- 
gular momentum is not imposed. 



3 COLLAPSE OF A ROTATING CLOUD 

INDUCED BY EXTERNAL COMPRESSION 

In this section we present the results of two simulations in- 
volving a core which initially has Po = 0.02. In the first 
simulation, the core is compressed slowly (4> = 3), and in 
the second it is compressed rapidly (<f> = 0.3). We limit the 
discussion here to a description of the density and velocity 
fields which develop on scales much larger than the central 
primary protostar or the rotationally supported disc which 
forms around it. A detailed discussion of the structure and 
evolution of the disc will be given in Section 2] Since the 
initial rotation energy is small (/3o = 0.02), rotation has lit- 
tle effect on the large-scale fields under discussion in this 
section, and most of the dynamical effects are the same as 
for the non-rotating cores analyzed in Paper I. 



3.1 Slow compression (<j> = 3, /3o = 0.02) 

In Fig. 0we show results for slow compression, = 3. Four 
times are shown: t = 0.560 Myr (thin full line) is well before 
the central primary protostar forms; t — 0.610 Myr (dot- 
ted line) is when the maximum density first reaches po, just 
before the central primary protostar forms; t = 0.611 Myr 
(dashed line) is after the central primary protostar has 
formed and the disc has just started to form; t — 0.615 Myr 
(dot-dash line) is about 4000 years after the disc starts 
to form. Plots (a) and (b) are log-log plots showing, re- 
spectively, the run of density along the equatorial plane 
(p(w,z — 0)) and the run of density along the polar axis 
(p{w = 0, z)); for reference, the thick full line on these plots 
shows the density of the singular isothermal sphere (here- 
after SIS), 



Psis — 



Of 



2-wGr 2 



(13) 



Plots (c), (d) and (e) are log- linear plots showing, respec- 
tively, the run of equatorial velocity (v m (w, z — 0)), the run 
of polar velocity (v z (w = 0,2)), and the run of azimuthal 
velocity in the equatorial plane (v$(w, z — 0)). 

In the outer parts of the core (r > 0.03 pc), the density 
profile is very close to the SIS, both along the equator, and 
along the pole. However, towards the centre, the equatorial 
density (p(w, z = 0)) gradually becomes larger than the SIS 
density, and the polar density {p(w = 0,z)) gradually be- 
comes smaller than the SIS density. The reasons for this are 
analyzed in the Appendix. 

In preparation for the analysis in the Appendix, Figure 
|5Ja) shows the equatorial density profile at t = 0.610, when 
the density first exceeds po, normalized to the density of the 
singular isothermal sphere (thin dashed line) . In addition, we 
have simulated the collapse of the same core (Bonnor-Ebert 
sphere with £ = 6), compressed at the same rate (0 = 3), 
but with no rotation (/3o = 0); the radial density profile 
obtained in this case, when the density first exceeds po at 
time t = 0.497, is shown as a dotted line. (The thick dashed 
line on this plot is defined in the Appendix.) The equatorial 
density profile of the non rotating core is higher than the 
density profile of the SIS, and the equatorial density profile 
of the rotating core is higher still, particularly towards the 
centre. 

On Figure the equatorial density increases abruptly 
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Figure 1. Slow compression (tj> 
with /3o = 0.02. Plot (a) shows log 10 [p(ui 
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Plot (b) shows log 10 [p(w = 0, z)/gcm — 3 ] 
Plot (c) shows v w (w,z = 0)/kms 



against logj^Mi/pc]. Plot (d) shows v z (w = 0, z)/kms" 
against log 10 [z/pc]. Plot (e) shows ti((i»,z = 0)/kms _1 against 
log 10 [w/pc]. Four times are shown: t = 0.560 Myr (thin full line) 
is well before the central primary protostar forms; t = 0.610 Myr 
(dotted line) is when the density first reaches po, j us t before the 
central primary protostar forms; t = 0.611 Myr (dashed line) is 
after the central primary protostar has formed and the disc has 
just started to form; t = 0.615 Myr (dot-dash line) is about 4000 
years after the disc starts to form. For reference, the thick full line 
on plots (a) and (b) shows the density of the singular isothermal 
sphere, C$/2nGr 2 (see Eq.[T31. 
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Figure 2. These plots show equatorial density profiles, p(w,z = 
0), normalized to the density profile of the singular isothermal 
sphere (see Eg. 1131 . for the two cases: (a) = 3, i.e. slow com- 
pression; and (b) cf> = 0.3, i.e. fast compression. The thin dashed 
line gives the equatorial density profile when the density first 
exceeds po, close to the moment the central primary protostar 
forms, and before the disc starts to form; in case (a) this mo- 
ment is t = 0.610 Myrs, and in case (b) it is 0.2585 Myr. The 
dotted line gives the equatorial density profile, again when the 
density first exceeds po close to the moment the central primary 
protostar forms, but for a non-rotating core (all the other initial 
conditions of the core are the same); in case (a) this moment is 
t = 0.497 Myr, in case (b) it is 0.225 Myr. The thick dahed line 
gives the product of the density of the non-rotating cloud (dotted 
line) and the factor 1 + (vg/C 3 ) 2 /2 (see Equation IA10I . 



(by a factor 10 to 20) inside r ~ 2x 10~ 4 pc at t = 0.611 Myr, 
and inside r~5x 10~ 4 pc at t = 0.615 Myr. This abrupt 
density increase marks the accretion shock at the outer edge 
of the growing disc. 

The inward velocity at the edge of the core is sim- 
ilar at the poles and around the equator, ranging from 
~ 0.07kms _1 at t = 0.560 Myr to ~ 0.11 kms" 1 at t > 
0.610 Myr. Towards the centre, the magnitude of the po- 
lar velocity \v z (w — 0,z)\ increases more rapidly than the 
magnitude of the equatorial velocity \v w (w,z — 0)|, due to 
centrifugal acceleration. Interior to 0.01 pc, \v z (w — 0, z)\ 
is approximately twice \v w (w,z = 0)|, until the material 
flowing inwards close to the equator encounters the outer 
boundary of the disc. At this point, \v w (w,z = 0) | decreases 
abruptly in the accretion shock at the disc boundary. The 
maximum value of \v w (w,z = 0)|, just before the material 
hits the disc boundary, is approximately constant at a value 
~ 0.85 kms -1 . In contrast, the maximum polar velocity in- 
creases constinuously. 

This simulation has been repeated with the procedure 
described in Sect. 12.61 which forces local conservation of an- 
gular momentum. The results are very similar, the only dif- 
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Figure 3. Same as Fig. Q but for fast compression (0 = 0.3) 
of a rotating cloud with j3o = 0.02. Four times are shown: 
t = 0.240 Myr (thin full line) is well before the primary pro- 
tostar forms; t = 0.2585 Myr (dotted line) is when the density 
first reaches poi just before the central primary protostar forms; 
t = 0.2593 Myr (dashed line) is after the central primary pro- 
tostar has formed and the disc has just started to form; and 
t = 0.2625 Myr (dot-dash line) is about 3000 years after the disc 
starts to form. The thick line on plots (a) and (b) is the density 
of a singular isothermal sphere. 



ference being that the maximum value of \v w (w, z = 0) | , just 
before the material hits the disc boundary, is ~ 0.75 kms - 
instead of ~ 0.85 km s -1 . 

3.2 Fast compression (0 = 0.3, = 0.02) 

In Fig. [3]we show the results for fast compression, = 0.3. 
Four time are shown: t = 0.240 Myrs (thin full line) is well 
before the central primary protostar forms; t = 0.2585 Myr 
(dotted line) is when the maximum density first reaches 
po just before the central primary protostar forms; t = 
0.2593 Myr (dashed line) is after the central primary pro- 
tostar has formed and the disc has just started to form; and 
t — 0.2625 Myr (dot-dash line) is about 3000 years after the 
disc starts to form. Some significant differences from the case 
(f> — 3 can be seen in the density and velocity profiles. 

First, due to the more rapidly increasing external pres- 
sure, the core is more compact; the cloud radius is about 
0.035 pc for = 3 whereas it is about 0.025 pc for = 0.3 
(see Figures Q and . Consequently, the equatorial density 
p(w, z = 0) is higher than for = 3 (typically by a factor ~ 
1.5-2). In contrast, the polar density is more or less the same 
as for the case = 3. From Fig. EI there appear to be two fac- 
tors contributing to the increase in equatorial density with 
decreasing 0. (a) The non-rotating cloud (dotted lines) has 
higher equatorial density for = 0.3 (middle panel) than for 
= 3 (upper panel), so the first factor depends only on the 
rate of compression, and not on the angular momentum, (b) 
The ratio of the rotating cloud density to the non-rotating 
cloud density is larger for = 0.3 than for = 3, so the 
second factor depends both on the rate of compression, and 
on the angular momentum. The reason for this is analysed 
in the Appendix. 

Second, the inward equatorial velocity \v w (w,z = 0)| in 
the outer parts of the core is greater for = 0.3 than for 
= 3. For example, the edge velocity is ~ 0.15 — 0.20 kms~ 
for = 0.3 (as compared with ~ 0.07 — 0.11 kms - for 
= 3); and at r = 0.01 pc, it is ~ 0.3 — 0.35 kms -1 for 
= 0.3 (as compared with ~ 0.18 — 0.20 kms -1 for = 3). 
The equatorial velocity profile is also flatter than for the case 
= 3, due to the compression wave. At t = 0.2625 Myr, 
there are large fluctuations in v w (w,z = 0) at small radii 
(r < 0.002 pc), due to the development of non axisymmetric 
modes. These will be discussed further in Section [I] 



4 FRAGMENTATION 

In this section we focus on what happens in the central parts 
of the cloud, and the details of the fragmentation process. 

4.1 Slow compression 

Fig. 0] illustrates the development of instability in the disc 
which forms around the central protostar, for the case = 3 
(slow compression) and f3o = 0.02. The lefthand column 
shows particle positions projected onto the z = plane. 
The righthand column shows logio [pi] plotted against equa- 
torial co-ordinate Wi, for each particle i having density above 
10 -14 gcm -3 , and the solid line shows p(w), the mean den- 
sity interior to radius w, as defined in Equation 1A3H . The 
central density first rises above po/3 at t = 0.6097 Myr, and 
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Figure 4. Disc instability in the central ~ 10 pc for = 3 
(slow compression) and /3o = 0.02. The lefthand column shows 
particle positions projected onto the z = plane. The righthand 
column shows log\o\pi\ plotted against equatorial co-ordinate 
Wi, for each particle i having density above 10 — 14 gem -3 , and 
the solid line shows p(u>), the mean density interior to radius 
w, as denned in Equation IA31I . Five timesteps are shown, 
t = 0.6117, 0.6127, 0.6137, 0.6147, 0.6177 Myr. The mass on top 
of the panels is in Mq. The density and velocity fields on larger 
scales are illustrated in Fig. No fragmentation occurs. 
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Figure 5. Disc instability in the central ~ 2 X 10 pc for i 



: 0.3 



(rapid compression) and /3o = 0.02. The lefthand column shows 
particle positions projected onto the z = plane. The righthand 
column shows logio[pi] plotted against equatorial co-ordinate 
Wi, for each particle i having density above 10 -14 gcm -3 , and 
the solid line shows p(u>), the mean density interior to radius 
w, as denned in Equation IA31I . Five timesteps are shown, 
t = 0.2599, 0.2608, 0.2619, 0.2628, 0.2646 Myr. The mass on top 
of the panels is in Mq. The density and velocity fields on larger 
scales are illustraated in Fig.|3] A second protostar condenses out 
of one of the spiral arms. 



the five timesteps shown correspond to 2000, 3000, 4000, 
5000, and 8000 years after this. The density and velocity 
fields on larger scales are illustrated in Fig. Q 

Panel (a) on Figure 2] shows the disk which forms in 
the centre of the core around the primary protostar. At this 
stage its central density is about 3 x 10 -12 gcm -3 and its 
edge density is about ten times smaller. It is bounded by an 
accretion shock, where the density falls by a further factor of 



ten. The disc has (5 ~ 0.36, but it is still apparently symmet- 
ric. The density throughout the disc has only just started 
to rise above p, and therefore it is only mildly unstable ac- 
cording to the analysis presented in the Appendix (Section 
IA4|l . Symmetry breaking is first evident a few hundred years 
later, by which time (5 ~ 0.40. 

In panels (b) and (c) of Figure 2] a strong two-armed 
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spiral pattern develops in the disc, due to spontaneous sym- 
metry breaking, but the arms do not sweep up sufficient 
mass to become gravitationally unstable, and they quickly 
wind up. This is very reminiscent of the numerical results 
reported by Durisen et al. (1986), and is due to the fact that 
the m = 2 modes are the first to become dynamically un- 
stable (Chandrasekhar 1969, Ostriker & Bodenheimer 1973). 
As in the numerical simulations of Durisen et al. (1986), the 
arms generate gravitational torques which transport angular 
momentum outwards through the disc, allowing the central 
parts to condense onto the primary protostar, and the outer 
parts to expand. However, the situation we simulate here dif- 
fers from that modelled by Durisen et al., because the discs 
in our simulations are accreting from an infalling envelope. 
This has two fundamental consequences. First, the mass at 
the outer edge of the disc is continuously replenished; this 
effect was described by Bonnell (1994) and by Whitworth 
et al. (1995), who showed numerically that this process will 
often lead to fragmentation. Second, the edge of the disc is 
bounded by an accretion shock which compresses the gas at 
the edge of the disk (see Fig. 0. 

We note that since the perturbations that lead to sym- 
metry breaking are numerical noise due to the initial particle 
distribution, the details of the structures that form - for ex- 
ample, the orientation of the spiral arms - are not identical 
for two differents realizations, i.e. two different initial par- 
ticle distributions representing the same macroscopic initial 
conditions. However, the statistical properties - such as the 
numbers, masses and orbits of protostars produced - do not 
significantly depend on the initial particle distribution, nor 
do they depend significantly on the numerical resolution. 

Panels (d) and (e) of Figure0|show the subsequent evo- 
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lution of the disc. The disc is replenished by infalling mate- 
rial, and a second strong two-armed spiral pattern develops, 
but again it fails to sweep up sufficient material to become 
gravitationally unstable. Instead, it generates gravitational 
torques which transport angular momentum outwards, al- 
lowing the inner material of the disc to accrete onto the 
central primary protostar and dispersing the outer material 
of the disc. No secondary protostar is formed. 



4.2 Fast compression 

Fig. |S] shows results for <f> — 0.3 (rapid compression) and 
/3o = 0.02. In this case the central density first rises above 
po/3 at t ~ 0.2578 Myr (see Figure |SJ, and the timesteps 
shown are 2000, 3000, 4000, 5000 and 7000 years after this. 
The principal effects of more rapid compression are (i) to 
drive material into the disc more rapidly, thereby building 
up the mass of the disc more quickly, and curtailing the time 
the disc has to stabilize itself by redistributing angular mo- 
mentum and accreting onto the central primary protostar; 
and (ii) to increase the density in the outer parts of the disc 1 . 
The result is that the central primary protostar is smaller 
and the disc fragments to produce a secondary protostar. 

This can be seen in Panel (a) of Figure which illus- 
trates the structure of the disc just before symmetry break- 
ing occurs. In comparison with the case cj> = 3, the disc here 
is both more massive and has a flatter density profile, i.e. the 
central density is lower and there is more mass in the outer 
parts of the disc. Consequently the density in the outer parts 
is significantly higher than the mean density, and the disc is 
gravitationally unstable according to Condition 1A33I . At 
this stage (3 ~ 0.48, and by the time symmetry breaking 
occurs, /3 ~ 0.51. 

The development of the instability is similar to the pre- 
vious case, but because the outer parts of the disc are denser, 
and the central primary protostar is less massive, the spiral 
arms are denser, more extended, and consequently more un- 
stable (than for = 3). Increasing the rate of compression 
does not simply accelerate the formation of the disc, but 
also reduces the time available for redistribution of angu- 
lar momentum by symmetry breaking, thereby generating a 
greater ratio of disc mass to primary protostar mass. 

Panel (d) of Figurc[^]shows the non-linear development 
of the spiral arms, and in Panel (e) a second object forms, lo- 
cated at x ~ —6 x 10~ 4 pc and y ~ 1 X 10~ 4 pc. At this stage 
(t = 0.2646 Myr), the mass of the central primary protostar 
is Mi ~ O.O8M0 (75% of which was already in the system at 
time t = 0.2599 Myr, i.e. before symmetry breaking started), 
and the mass of the newly-formed secondary is ~ O.O2M0 
(of which 40% was in the system before symmetry breaking 
started). 

The systematic differences between the cases (j> — 3 
(slow compression) and (f> = 0.3 (fast compression) are fur- 
ther illustrated in Figure|S] which shows the mass of gas with 



We stress that this latter effect is not because the Mach num- 
ber of the accretion shock is higher - the speed with which ma- 
terial flows into the accretion shock at the edge of the disc is not 
strongly dependent on the rate of compression - but because the 
density of the material flowing into the shock is higher for faster 
compression. 
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density in excess of various representative thresholds (10po, 
3po, Po, Po/3, po/10 and p/100), as a function of time. Sym- 
metry breaking occurs at t = 0.612 Myr for <j> = 3, and at 
t — 0.2632 Myr for (f> = 0.3. The maxima exhibited by the 
curves for intermediate thresholds (po, thin full line; po/3, 
thick full line) are due to expansion of the disc, caused by 
symmetry breaking and redistribution of angular momen- 
tum. Panel (a) of Figure |S| shows that for <f> — 3, when sym- 
metry breaking occurs, the density in the disc is relatively 
low, and a large fraction of the system mass is already in 
the central primary protostar, thereby stabilizing the disc. 
Conversely, Panel (b) shows that for (f> = 0.3, when sym- 
metry breaking occurs, the density in the disc is relatively 
high, and a much smaller fraction of the system mass is in 
the central primary protostar. Panel (b) also shows how the 
transport of angular momentum in the disc is accelerated 
once the secondary protostar forms at t — 0.2632 Myr. 

4.3 Faster compression: ring formation 

Figure □ shows results for 4> = 0.1, A) = 0.02. The be- 
haviour is very different from the previous cases. Even be- 
fore the maximum density approaches po = 10~ 13 gem -3 , a 
ring forms, and there is no central primary protostar. Ring 
formation is attributable to a combination of factors. 

First, there is a centrifugal barrier, and this creates the 
rarefaction at the centre of the ring. The dynamics of ring 
formation due to a centrifugal barrier have been analysed 
by Tohline (1980), on the basis of pressureless collapse in an 
external potential. Bonnell & Bate (1994) have also noted 
the transition from disc formation to ring formation, as the 
speed of collapse increases. In their case the speed of col- 
lapse was increased by reducing the initial ratio of thermal to 
gravitational energy in the cloud. Cha & Whitworth (2003) 
have explored the influence of differential rotation on ring 
formation. In order to understand better the origin of the 
ring, we have monitored where the material impinging on the 
centre of a core originates, and we find the following distinc- 
tion. For relatively slow compression (c^>0.3, sections 4.1 
& 4.2), the material which first impinges on the centre of 
the core was initially concentrated near the rotation axis (z- 
axis). Therefore this material has very low specific angular 
momentum (as compared with material originating further 
from the rotation axis). This is why it reaches the centre 
first (it experiences least centrifugal acceleration than ma- 
terial originating further from the rotation axis), and why on 
reaching the centre it can stay there to form the central pri- 
mary protostar. In contrast, for faster compression (<^><0.1, 
this section), the compression wave is stronger, and drives 
material into the centre more isotropically. As a result, most 
of the material impinging on the centre originates far from 
the rotation axis and therefore has too much angular mo- 
mentum to reach the centre, so there is a central rarefaction 
- and hence a ring is formed. 

Second, the material delivered into the outer parts of 
the nascent disc is compressed to high density by the accre- 
tion shock at the edge of the disc. In order to confirm this 
effect, we have extracted the ring displayed in panel (b) of 
Figure 7 (i.e. we have selected particles having density larger 
than po/10) and we have then let the ring evolve in isolation 
whilst enforcing axisymmetry. The ring quickly settles into 
an equilibrium in which ~ 17% of its mass has density below 
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Figure 7. Disc instability in the central ~ 10 pc for <j> = 0.1 
(very rapid compression) and /3rj = 0.02. The lefthand col- 
umn shows particle positions projected onto the z = plane. 
The righthand column shows log\o[pi\ plotted against equato- 
rial co-ordinate Wi, for each particle i having density above 
10~ 14 gcm — 3 , and the solid line shows p(ui), the mean density 
interior to radius to, as defined in Equation IA31I . Five timcstcps 
are shown, t = 0.1688, 0.1704, 0.1712, 0.1717, 0.1744 Myr. The 
mass on top of the panels is in Mq. Fragmentation occurs via 
a ring which initially breaks up into three pieces; two of these 
merge, and the end-result is two protostars. 



po/10, and this mass carries ~ 33% of the angular momen- 
tum of the ring. We conclude that it is the accretion shock 
which gives the ring a sharply defined outer boundary, and 
maintains its high mean density and high specific angular 
momentum. 

Third, once the ring density becomes higher than the 
mean density, its gravity starts to attract infalling material 
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away from the centre, thereby further enhancing its den- 
sity contrast (see Panel (b) of Figure Q t = 0.1704 Myr) 
The ring is established so quickly that there is insufficient 
time for the symmetry-breaking instabilities which could re- 
distribute angular momentum and deliver material into s 
central primary protostar. 

Self-gravitating rings are very unstable to non- 
axisymmetric instabilities (Ostriker 1964; Norman & Wil- 
son 1978), and within a few orbital periods of its forma- 
tion it breaks up into three massive fragments (Panel (c) 
of Figure t = 0.1712 Myr). At the same time these frag- 
ments start to interact dynamically (Panel (d) of Figure |7J 
t = 0.1717 Myr). Two of them merge, and the end result 
is a binary (Panel (e) of Figure |7| t = 0.1744 Myr). The 
object located at x = —0.1 x 10~ 3 pc, y = —0.3 x 10~ 3 pc 
on Panel (e) of Figure Q has mass ~ O.O8M0; the object at 
x = +0.3 x 10~ 3 pc, y = +0.6 x 10~ 3 pc has mass ~ O.O6M 
Both objects are still accreting and have small discs with 
spiral patterns around them. Ring fragmentation tends tc 
produce objects of comparable mass (as here), in contrast 
to disc instability where the secondary protostars formed foi 
4> — 0.3 tend to be four or five times less massive than the 
primary. 

In order to demonstrate the importance of the ram pres- 
sure of the accretion shock, we have performed a simple 
numerical experiment. We first extract the ring structure 
displayed in panel (b) of Figure |7J (i.e. we select particles 
having a density larger than po/W). Then we let this ring 
evolve in isolation, first with no external pressure, and sec- 
ond with an external pressure equal to the average therma 
pressure of its particles (i.e. comparable to the ram pressure 
of the infalling gas). In the first case (no external pressure 
lefthand column of Figure [8J, there is no permanent frag- 
mentation. Transient structures develop in the ring, but, be- 
cause they lack a confining pressure, they are diffuse, anc 
they merge to form a single central protostar. The rest oi 
the material ends up in an expanding disc. A similar result 
has been reported by Bonnell (1994), who finds that the 
disc which forms around his central primary protostar only 
fragments if envelope material continues to fall in onto the 
disc. In the second case (external pressure approximately 
equal to ram pressure, righthand column of Figure |HJ , the 
evolution is broadly similar to that presented in Figure |7| 
The ring breaks up into three fragments, and subsequently 
two of them merge. This indicates that the confinement of 
the ring by the ram pressure of the infalling gas plays an 
important role in promoting fragmentation. 



4.4 Stronger rotation 

In order to investigate the effect of higher rotation on frag- 
mentation, we have performed numerical simulations with 
/3o = 0.05. The results for = 3 and /3o = 0.05 are pre- 
sented in Figure |U] In this case the central density first rises 
above po/3 at t — 0.7237 Myr, and the timesteps shown are 
2000, 3000, 4000, 5000 and 7000 years after this. The prin- 
cipal effect of higher angular momentum is to increase the 
mass and extent of the disc at the expense of the central pri- 
mary protostar (relative to the case cj> = 3, /3o = 0.02), so the 
disc is more unstable. As a result, the spiral arms become 
self-gravitating and condense out to produce two secondary 
companions. 
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Figure 8. Evolution of the ring structure displayed in panel (b) 
of Figure when the envelope is removed. The lefthand column 
shows the evolution when there is no external pressure. With 
no external pressure, a single central protostar forms, the outer 
parts of the ring expand, and there is no fragmentation. The 
righthand column shows the evolution when the external pressure 
is approximately equal to the ram pressure of the infalling gas. 
In this case, the evolution is very similar to the full simulation 
displayed in Figure 171 



We infer that, in the parameter space that we have 
explored here (and that we expect to be representative of 
real star- forming cores), both rotation (higher /3o) and rapid 
compression (higher (j>) promote fragmentation. This is in 
accordance with the analysis in the appendix. For higher <j> 
and Po the material impinging on the accretion shock at the 
edge of the disc has higher density, and therefore the den- 
sity in the outer parts of the disc is higher. In addition, for 
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5 SUMMARY AND DISCUSSION 
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Figure 9. Disc instability in the central ~ 10 pc for = 3 
(slow compression) and 0o = 0.05 (above average rotation). The 
lefthand column shows particle positions projected onto the z = 
plane. The righthand column shows logio[pi] plotted against 
equatorial co-ordinate u>i , for each particle i having density above 
10 _14 gcm — 3 , and the solid line shows p(u>), the mean density 
interior to radius w , as defined in Equation IA31I . Five timesteps 
are shown, t = 0.7256, 0.7268, 0.7277, 0.7287, 0.7307 Myr. The 
mass on top of the panels is in Mq. Fragmentation of the disc 
arouund the primary protostar produces two further protostars. 



higher (j> and /3q , this material is delivered to the outer parts 
of the disc more rapidly, so there is less time for disc ma- 
terial to redistribute angular momentum, and the mass of 
the central primary protostar is lower. Therefore, as demon- 
strated in the appendix, the accretion and local time-scales 
are reduced relative to the global time-scale, and the disc is 
more unstable against fragmentation. 



We have investigated the effect of compression on the col- 
lapse and fragmentation of a rotating core. 

Most of the conclusions of Paper I, concerning the large- 
scale dynamics of the collapsing core, are still valid. In par- 
ticular, the increase in external pressure drives a compres- 
sion wave into the core. The compression wave leaves in 
its wake a velocity field very similar to those recently in- 
ferred for prestellar cores, from observations of asymmetric 
molecular-line profiles. Tafalla et al. (1998) estimated inflow 
at ~ 0.10 km s" 1 in the outer layers of L1544, and Williams 
et al. (1999) found inflow at ~ 0.08 km s -1 further in. Lee 
et al. (1999) detected inflow velocities ranging from 0.04 
to O.lOkms -1 in several other prestellar cores. Our models 
generate comparable inflow velocities for <j> in the range 10 
to 1, i.e slow to intermediate compression. Interestingly, in 
regions like Perseus and p Ophiucus, where star formation is 
triggered more violently, cores with significantly higher infall 
velocities have been observed. For example, Di Francesco et 
al. (2001) detected velocities of ~ 0.5 to 0.7kms _1 in the 
Class protostars of NGC 1333 IRAS 4, suggesting that 
very fast compression may have occurred in this region. 

The fundamental difference from Paper I is that because 
the cores we simulate here have rotation, the inflowing ma- 
terial does not all converge directly onto the central primary 
protostar (CPP). Instead, a large fraction of it first collects 
in an accretion disc around the CPP, and instabilities in 
the outer parts of this disc may then lead to fragmentation 
producing additional protostars. As noted by Larson (2002) , 
the critical factors determining the stability of the outer disc 
are (i) its mass and density, and (ii) the speed with which it 
is assembled: a quickly assembled, massive, dense outer disc 
is unstable against fragmentation. 

We show numerically - and, in the Appendix, semi- 
analytically - that the density in the outer disc is larger 
for faster compression, and for larger /3q. These effects may 
already have been observed. In cores belonging to the rel- 
atively quiescent star formation region Taurus, the inflow 
velocities are compatible with slow compression, and the 
density is close to the density of the singular isothermal 
sphere (SIS). Conversely, in cores belonging to more active 
star formations regions, the inflow velocities appear to be 
supersonic, and the densities are about one order of mag- 
nitude higher than the density of the SIS (Motte & Andre 
2001; Andre et al. 2003). This accords with the predictions 
of the analysis in the appendix. 

Except in the case of very rapid compression, the low 
angular momentum material arriving in the centre of the 
core goes directly into the CPP, and the high angular mo- 
mentum material forms an accretion disc round the CPP. If 
the compression is more rapid, this has a number of effects 
which tend to render the disc more unstable. First, the ma- 
terial accreting onto the outer parts of the disc arrives with 
higher density, and therefore the density in the outer disc 
- following compression in the accretion shock at the edge 
of the disc - is also higher. Second, material is delivered 
into the outer parts of the disc more rapidly, and there- 
fore these outer parts become more massive. Third, there 
is less time for the gravitational torques associated with 
symmetry-breaking instabilities in the disc (the same insta- 
bilities which lead to fragmentation) to redistribute angular 
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momentum and thereby facilitate the continuing growth of 
the CPP. The combination of a massive, dense outer disc 
and a low-mass CPP makes the outer disc unstable against 
fragmentation, spawning secondary protostars with masses 
typically four or five times lower than the CPP. 

For very rapid compression there is no CPP; all the ma- 
terial flows into the disc, and it is so concentrated towards 
the edge that it is more accurately described as a ring. The 
ring then fragments into two or three protostars of compa- 
rable mass. 

For more rapid rotation (f3o = 0.05) the outer disc is 
even more massive in comparison to the CPP, even more 
extended, and therefore even more prone to fragment. 

In the appendix we analyse the structure of the in- 
flowing envelope, and its consequences for the stability of 
the central disc. This analysis explains why higher rotation 
(large /3o) and more rapid compression (small d>) promote 
fragmentation and the formation of multiple protostars. 
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In the appendix we analyse the structure of the in- 
flowing envelope, and its consequences for the stability of 
the central disc. This analysis explains why higher rotation 
(large /3q) and more rapid compression (small cj>) promote 
fragmentation and the formation of multiple protostars. 



APPENDIX A: ANALYTICAL PREDICTIONS 

In this appendix, we develop an analytic description of the 
structure of the inflowing envelope and its consequences for 
the stability of the disc which forms around the central pri- 
mary protostar (or, in the case of very rapid compression, 
the ring structure which forms around a central rarefaction), 
with a view to understanding how disc fragmentation is af- 
fected by changes in the initial rotation (/3o) and the rate of 
compression (0) . The formation and evolution of a disk em- 
bedded in a rotating and collapsing core has already been 
investigated and elegantly described in seminal papers by 
Cassen & Mossman (1981) and Stahler et al. (1994). How- 
ever, these authors assume that the core collapses according 
to the inside-out model of Shu (1977), starting from a sin- 
gular isothermal sphere. Here we wish to consider the case 
of dynamically triggered collapse, which is necessarily from 
the outside-in. 

Four features of the collapse are crucial to this discus- 
sion. First, the polar density profile p(w = 0, z) is close to the 
density profile of the SIS, whereas the equatorial density pro- 
file p(w,z — 0) is significantly higher. Second, the maximum 
value of the inward equatorial velocity \v w (w, z = 0)| ma x in- 
creases monotonically with time until the disc forms, af- 
ter which it is approximately constant. Third, this asymp- 
totic maximum equatorial velocity (which determines the 
strength of the accretion shock at the edge of the disc) de- 
pends only weakly on the rate of compression (j>. Fourth, the 
strength of the accretion shock at the edge of the disc has 
an important influence on the stability of the disc. By ana- 
lyzing these effects, we can estimate the different timescales 
controlling fragmentation, and hence interpret the results 
reported in Section 2] 



Al Density profile 

According to the numerical results displayed on Fig. [5] the 
equatorial density profile p(w, z = 0) depends on the rate of 
compression (f> and on the initial rotation j3o. 



A 1.1 The effect of external pressure 

The first effect (rapid compression leading to large equato- 
rial density) can be understood qualitatively by reference 
to the self-similar solutions studied by Whitworth & Sum- 
mers (1985). In these solutions, the density at large radius 
converges to Uoo/r 2 . Whitworth & Summers (1985) show 
that the stronger the compression wave being driven into the 



core, the faster the collapse and the higher H ro . The slow- 
est collapse corresponds to Shu's inside-out collapse from 
a singular isothermal sphere (Shu 1977) and has Uoo — 1. 
The Larson-Penston solution (Larson 1969, Penston 1969) 
corresponds to collapse from a centrally flat density profile, 
induced by a strong compression wave, and has ~ 7. 

We therefore assume that the density profile can be ap- 
proximated by 

P(r)*±, (Al) 

where A is a constant. This assumption is justified both 
by the numerical results presented in Figure [5] and by the 
asymptotic form of the similarity solutions obtained by 
Whitworth & Summers (1985). Significant departures from 
p(r) oc r~ 2 are confined to the inner parts of the core which 
contain very little of the total mass. 

If R c is the core radius, then the core pressure at this 
point must be equal to the external pressure, P ex t, i.e. 



Mass conservation requires 
4nR c A = M C - M,, 



(A2) 



(A3) 



where M c is the initial core mass and M* is the mass of the 
central primary protostar. As long as M» <C M c , as it is 
when the disc first forms, M* can be neglected, so 

2/3 



.4 



1/3 / 



M c 



(A4) 



Recalling that the density of a singular isothermal sphere is 



Psis = 



2wGr 2 
we can write 
p(r) 2irGA 



1/3 



(A5) 



(A6) 



Psis(r) 

Here P ~ Cf/G 3 M C 2 is the pressure at the boundary of the 
core before compression starts. From the numerical results 



lO- 



gcm 3 , giving p/psis ^ 1-42; whilst 



gem 3 , giving p/p S rs ^ 1-93. The 



for (j> = 3, p cxt 
for (f> = 0.3, p cxt — 10 
dotted lines on Figure [5] demonstrate that these predictions 
are in good agreement with the numerical results. 

A 1.2 The effect of rotation 

The similarity equations describing self-gravitating collapse 
have been extended to include rotation by Saigo & Hanawa 
(1998) for disk geometry, and by Hennebelle (2003) for fila- 
mentary geometry. However, the geometry assumed in these 
treatments is significantly different from the geometry that 
we are considering here. In order to investigate analytically 
the effect of rotation on the density profile, we simply as- 
sume that the core is close to equilibrium, so that in the 
equatorial direction 



c 2 dp 



p dw w dw 



(A7) 



If we now neglect departures from spherical symmetry, and 
substitute 



p(w,z = G)~—. 



(A8) 
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so that d<&/dw ~ —iirGA/w, we obtain 



.4 



2txG 



1 + 



2CI 



and hence 
p(w,z = 0) 
psis(w) 



2d 



(A9) 



(A10) 



If vg = 0, the SIS density profile is recovered. In general, ve 
is finite, but not constant, and so A is not constant either. 
However, the variation of Ve with w is always much weaker 
than w~ 2 , so to a first approximation we can treat A as con- 
stant. Equation (IA10I explains why the equatorial density 
profile p(w, 2 = 0) is higher than the SIS density profile, and 
why the increase is greatest in the inner parts of the core, 
where ve is greatest (see Figure 0. 

In order to test the predictions of Equation IA10I , Fig- 
ure |5| compares the density profiles obtained in simulations 
with <f> = 3 or 0.3 and /3o = 0.02 (thin dashed lines) with 
the density profiles obtained in simulations with no rotation 
(dotted lines), but then multiplied by the factor 1 + Vg/2C 2 
from Equation l|A10(l to give the thick dashed lines. The 
comparison is made at the time when the maximum density 
first reaches po, and values of ve are taken from the rotating 
simulations represented by the thin dashed lines. In general 
the agreement between the two dashed lines (thin and thick) 
on Figure |21 is good, particularly for large values of (j>. 

Combining Equations l|A6^ and I|A1Q0 , we propose that 
by the time the central primary protostar forms, the equa- 
torial density profile can be approximated by 

p(w,z = 0) = 5 Psis (w) = -^J^, (All) 



2-kGw 2 



(A12) 



S is the factor by which the density in the inflowing gas 
in the envelope is enhanced by the combined effect of com- 
pression and rotation. In the outer parts of the envelope, the 
rotational velocity is normally very small compared with the 
sound speed, and so there the overdensity (relative to the 
SIS) is dominated by the effect of compression. However, in 
the inner parts of the envelope the rotation velocity becomes 
large, compared with the sound speed, and both effects are 
then important. 

A2 Infall velocity 

Now consider a parcel of gas, falling inwards onto the disc. 
Its acceleration is the sum of gravitational and centrifugal 
terms; we neglect the effect of thermal pressure. We assume 
that the mass of the disk contained within radius w is given 
by 



TTp W £, 



(A13) 



where po is the initial density of the core (assumed to be 
uniform, for simplicity), wo is the initial position of the gas 
parcel which is now (at time t) located at w, and £ is the 
initial height of the cylinder that has now flattened into the 
disc out to radius w. Let flo be the initial angular speed of 
the core. Then the equation of motion for the parcel is 

d 2 w 



FGM U 



+ 



(A14) 



where V is a geometrical factor of order unity. Equation 
can be integrated to give 

fdw\ 2 2hYG P q£w 2 wfa 2 , a 

-77 5-+ V 0, ( A15 ) 

V at J w w z 

where vo is the initial velocity and depends on wo and (5. 
It follows that the maximum equatorial velocity, 



\v w (w,z = 0)\ B 



is reached at 



«W*-^, (A16) 

and this is in effect the position of the accretion shock at 
the edge of the disc, where the parcel of gas that we are 
following accretes onto the disc. Therefore the edge of the 
disc is at 



Wedge ^ • (A17) 

7rl Gpot 

Combining Equations HA15II and HA16L the maximum 
inward equatorial velocity reached by the parcel of gas as it 
impinges on the accretion shock at the edge of the disc is 

dw 



dt 



fTnGpoly a 



n 

TitGpgt 
too 



1/2 



(A18) 



Equation <A18> shows that t> a cc depends only on Eq = 
pol (i.e. the initial cloud surface density or some fraction 
thereof) and not on wo- By the time the disk forms, the 
cloud is significantly flattened and this quantity is almost 
constant. Consequently, the variation of « aC c is expected to 
be small. 

Substituting from Equation 1 A 1711 in Equation 1A13II . 
we obtain an expression for the total mass of the disc 



Afdisc 



n 2 

"0 



(A19) 



We reiterate that these equations are valid only if the 
thermal pressure can be neglected. In particular, if f3a is 
too small then the gas parcel becomes adiabatic (p > po) 
before it reaches the disc. In the cases treated here, the gas is 
still isothermal when it first encounters the accretion shock 
at the edge of the disc (since its density is ~ O.lpo _ see 
Figures and |SJ . Thus neglecting the thermal pressure is 
acceptable, and the ram pressure of the gas flowing into the 
disc (~ pv 2 cc ) is larger than its thermal pressure (~ pC'o)- 

Since the values of w acc obtained in Section 3 appear to 
depend only weakly on the rate of compression (ij>), we infer 
that u acc 3> vo- Therefore neglect of vo in the final form of 
Equation IA18I is justified. It follows that for two different 
values of /3o, 

1/2 

(A20) 



IV 



(since Qo oc P^ 2 )- 

In order to compare the predictions of Equation 1A18I 
with the numerical results, we need to estimate the combi- 
nation pol/£lo- For the uniform-density, uniformly rotating 
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cloud from which our initial conditions are generated, we 
can write 



pol_ < 3M e 



(A21) 



where M c is the core mass and R c the core radius. Combin- 
ing Equations <A18t and 1A21II . we obtain 



3rGA/ c 



1.2kms 



(A22) 



Although this inequality has been derived assuming a 
uniform-density, uniformly rotating core, it is still valid for 
our simulations. That is because the stretching which cre- 
ates the initial conditions for our simulations (by converting 
a uniformly-rotating uniform-density core into a differen- 
tially rotating BE sphere) conserves both the core mass M c , 
and the specific angular momentum R 2 Qo- From Figures Q 
and|3]we see that in the simulations v acc ~ 0.8 ±0.1 kms -1 . 
We note that Equation IA22I gives an upper limit on w acc 
because £<R C , and because we have neglected the thermal 
pressure in deriving Equation <A18t . 

In order to compare the predictions of Equation <A20t 
with the numerical results, we have performed a simulation 
with /3o = 0.05 and cf> — 3 in which we find v aC c(0.05) ~ 
0.55kms _1 s, as compared with v acc (0.02) ~ 0.85 km s -1 in 
the simulation with /3q = 0.02 and </> — 3, giving a ratio 



« acc (0.05) 



; (0.02) 



0.65. 



simulation 



If we can neglect variations in T and t, Equation HA20t 
diets a ratio 



(A23) 



pre- 



=(0.05) 



=(0.02) 



/0.02\ 

Vao5/ 



1/2 



0.63 . (A24) 

analysis 

In view of all the approximations and assumptions made in 
deriving this result, the extreme closeness of the agreement 
between Eqn. I|A23(I and Eqn. l|A24fl must be somewhat for- 
tuitous, but it suggests that our analysis is a reliable guide 
to trends. 

Finally we can show that the tangential velocity of the 
parcel of gas which is about to impinge on the edge of the 
disc, iitang = vg(w — w c d S a, z = 0) should be approximately 
constant. The specific angular momemtum of the parcel is 
WqQo, so its tangential velocity is 

w 2 f3 nFGp £ 



'f'ta 



™cdg 



(A25) 



where we have obtained the second expression on the right- 
hand side of Equation 1A25I by substituting from Equation 
<A17t . Comparing Equation JA25t with Equation (IA18II , we 
see that 



i'ta 



(A26) 



and hence «t ang is approximately constant like fJ aC c- This is 
confirmed by the numerical results. For <f> — 3, the sim- 
ulations give Wedge = 2 x 10~ 4 pc, Wtang = 0.95 kms -1 
and Vacc = 0.87kms _1 at t = 0.611 Myr; and Wedge = 
5 x 10~ 4 pc, utang = 1.00 kms -1 , and u aC c = 0.85 kms -1 
at t — 0.615 Myr. For (j> = 0.3, the simulations give u> e dg C = 
3 x 10 _4 pc, Vt&ag = 1.00 kms -1 and u a cc = 0.69 kms -1 
at t — 0.259 Myr; and w c d gc = 10 X 10 -4 pc, utang = 
0.90 kms" 1 , -Oacc = 0.70 kms -1 at t = 0.262 Myr. 



A3 Accretion shock 

In order to analyze the accretion shock at the boundary 
of the disc, we define p e d gc to be the density just inside 
the edge of the disc (i.e. the post-accretion-shock density, 
p(w — Wedge — e, z = 0), where 2e is the shock thickness), 
Pace to be the density just outside the edge of the disc (i.e. 
the pre-accretion-shock density, p(w — w c dgc + e, z = 0)), and 
Vshock to be the outward equatorial velocity of the shock 
relative to the centre of the core. u aC c is the inward equato- 
rial velocity of the gas impinging on the shock at the edge 
of the disc (see Equation IA18II . Mass conservation requires 
Pcdgc^shock ^ Pace iiacc, and hence v s hock <S "an- Thus the 
velocity of the infalling gas in the shock frame is ~ v acc , and 
as long as the shock can be treated as isothermal, we can 
write 



Ped 



Pace 



\ Cs J ' 

From Equation IIA12I) . we have 
SC* 



Pace 



2^Gw c 2 dgc 



and so 



Pcdgc — 



5v 



ic.mnx 



2-kGwI 



(A27) 



(A28) 



(A29) 



(In the simulations presented in this paper, v acc /C a ~ 4 and 
so pedge should be ~ 16p a cc- This is corroborated by Figures 
IBI3liand0) 

A4 Time scales 

Disc fragmentation is an extremely non-linear process, gov- 
erned both by the intrinsic structure and evolution of the 
disc, and its interaction with the infalling material. Given 
the complexity of this interaction, the formulation of a pre- 
cise analytic criterion for fragmentation is probably impos- 
sible. However, useful insights can be gained by evaluating 
and comparing the timescales for competing processes. 

The global gravitational time scale for the disk is related 
to the orbital angular frequency, 

(GM W \' 1/2 f4nGp(w)^ ' " 

tglobal * { — ) * 

where we have introduced 
3M m 



(A30) 



(A31) 



P{W) 4^3 
the mean density interior to radius w. 

Similarly, the local gravitational time scale of the disk 
is related to the local Jeans frequency, 

ilocai = (4ttGpW)- 1/2 , (A32) 

where p(w) is the local density in the disc. 

The condition for instability is then that the local 
gravitational timescale be less than the global gravitational 
timescale, or equivalent ly 

\ 2 -/ \ 

tlocal \ PW 



^global 



p(w] 



< 1. 



(A33) 



(We note that this is essentially the same as Toomre's crite- 
rion, both for Keplerian discs, and for self-gravitating discs.) 
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To estimate Condition HA33t at the edge of the disc, 
we put p(w) — > pcdgc using Equation 1A29II and p(w) — > 
p(wedge) — 3Mdisc/47T , u;g d g C . Condition l|A33(l then becomes 



age 

pedge ^^cdge^acc 



p(Wedgc) 3GM dis 



> 1 . (A34) 



2TS _ 2T 

"3" ~~ T 



Finally, substituting for Mdisc, w e dge and v aC c from Equa- 
tions 1A19L \ A 1711 and 1A18L the condition for instability 
1A34I becomes 

(^)'"( 1 + |j)>i, 

where we have obtained the last expression by substituting 
for 8 from Equation 1A121 . This form of the condition for 
instability explains why more rapid compression (smaller 
4>) and more rapid initial rotation (larger /3q) both make 
the disc more unstable against fragmentation, by delivering 
higher density at the edge of the disc, i.e. higher 8. 

Another important time scale is the accretion 
time scale, taccrction = M d i sc /Mdisc- Putting M disc = 
27rw; e dgcftpacc^acc, where h is the vertical thickness of the 
layer of material flowing into the edge of the disc and PaccVacc 
is the flux of matter into the disc, we obtain 

^accretion — ~ i ■ 

(A36) 

27rW ed ge'iPacc«acc 

Then substituting for Mdi sc , Wedge , Pace and « acc from Equa- 
tions llATilll . llATFll . jHHt and llA"l"8t . Equation fKM re- 
duces to 

f-accrction t^i /)T-y^i9? " i } 

iaccrction is the timescale on which mass and angular mo- 
mentum are added to the disc, and it should be compared 
with t g i bai which is the minimum timescale on which mass 
and angular momentum can be redistributed within the disc. 
Substituting for tgiobai from Equation IA30I . and again us- 
ing Equations 1A13II and 1 A 1611 to eliminate Af d isc and w e d g e, 
we obtain 

^accretion _ TvCppilUQ . , 

tgiobai ~ rv*hCS8 ' [ ' 

A small value for this ratio implies an unstable disc, l/h is a 
geometrical factor, related to the cloud flattening, and is not 
easily calculated. Setting this factor aside, Equation l|A38fl 
implies firstly that large 8 (i.e. rapid compression and/or 
rapid initial rotation) promotes fragmentation, and secondly 
that small wo promotes fragmentation (i.e. fragmentation is 
more likely during the early stages of disc formation). 



